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Abstract 


A  three  dimensional  hydrostatic  finite  volume  ocean  model  is  developed  solving  the  integral 
dynamical  equations.  Since  the  basic  (integral)  equations  are  solved  for  finite  volumes  rather  than 
grid  points,  the  flux  conservation  is  easily  enforced  even  on  arbitrarily  meshes.  Both  upwind  and 
high-order  combine  compact  schemes  can  be  incorporated  into  the  model  to  increase  computational 
stability  and  accuracy.  For  abrupt  topography,  a  terrain-following  grid  discretization  is  designed  to 
reduce  computational  errors  such  that  the  four  lateral  boundaries  of  each  finite  volume  are 
perpendicular  to  x  and  y  axes,  and  the  two  vertical  boundaries  are  not  purely  horizontal.  This  grid 
system  reveals  a  superior  feature  than  z-  and  a  -coordinate  systems.  The  accuracy  of  this  model 
was  tested  in  this  study. 

Key  Words:  Finite  Volume,  Crystal  Grid,  Integral  Equations,  Hydrostatic  Balance,  Terrain- 
Following  Coastal  Model. 
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1.  Introduction 


Four  different  schemes  are  available  to  solve  partial  differential  equations  numerically:  (1) 
spectral  or  spectral  transform,  (2)  finite  difference,  (3)  finite  element,  and  (4)  finite  volume. 
Because  of  the  lateral  boundaries  with  complicated  shape,  the  ocean  basin  is  inherently  ill  suited  to 
the  spectral  technique.  The  finite  element  method  has  been  applied  to  2D  barotropic  problems 
extensively  such  as  tides  and  storm  surge  (Foreman  et  ah,  1993;  Le  Provost  et  al.  1994)  due  to  the 
flexibility  in  adapting  the  grid  locally  to  any  desired  resolution,  but  hasn’t  been  applied  to  3D 
baroclinic  problems  only  until  recent  years  (e.g.,  Lynch  et  al.,  1996).  A  principal  problem  of  this 
method  appears  to  be  the  mass  conservation.  While  globally  this  conservation  is  assured,  it  may  not 
conserve  the  mass  locally  (Kantha  and  Clayson,  2000). 

The  finite  difference  method  that  transforms  the  partial  differential  equations  into  difference 
equations  at  grid  points  is  commonly  used  in  regional  and  basin-scale  ocean  modeling.  Let  ( x ,  y) 
and  z  represent  the  horizontal  and  vertical  directions.  Various  finite  difference  models  use  different 
vertical  coordinates  such  as  z-coordinate  (e.g.,  Bryan,  1969),  terrain  following  cr  -  and  s- 
coordinates  (e.g.,  Blumberg  and  Mellor,  1987;  Song  and  Haidvogel,  1994),  and  isopycnal 
coordinate  (Bleck  et  al.,  1992).  The  solutions  of  the  finite-difference  models  are  valid  only  at  the 
grid  points.  For  coastal  oceans,  the  finite-difference  models  usually  use  the  terrain  following  cj - 
coordinate  and  have  large  truncation  errors  at  steep  topography  that  is  caused  by  horizontal 
pressure  gradient  errors.  Much  work  has  been  conducted  to  improve  the  accuracy  of  the  a  - 
coordinate  finite-difference  models  (e.g.,  Gary,  1973;  Beckman  and  Haidvogel,  1993;  Mellor  et  al., 
1994;  Chu  and  Fan,  1997,  1998,  1999,  2000,  2001). 

The  finite  volume  method  that  transforms  the  partial  differential  equations  into  integral 
equations  at  finite  volumes  has  yet  been  popular  in  ocean  circulation  modeling  and  simulation. 
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However,  the  conservation  is  easily  enforced  even  on  arbitrary  grids  because  the  integral  equations 
link  the  temporal  variability  of  the  dependent  variables  for  the  volume  to  the  fluxes  across  the 
boundary  of  that  volume  (Kobayashi,  1999;  Henneline,  2000).  This  leads  to  the  volume  setup  very 
flexible  that  makes  the  finite  volume  method  invaluable  especially  in  the  abrupt  topography. 

In  this  paper,  we  present  the  formulation  and  preliminary  test  of  the  finite  volume  ocean 
model  (VOM).  The  outline  of  this  part  is  as  follows:  A  description  of  the  dynamic  and 
thermodynamic  integral  equations  is  given  in  section  2.  A  depiction  of  the  finite  volume 
discretization  with  crystal  grids,  flux  computation,  and  explicit  finite  volume  scheme  is  given  in 
sections  3.  The  preliminary  model  test  case  is  depicted  in  section  4.  The  comparison  between  the 
finite  difference  and  finite  volume  methods  is  discussed  in  section  5.  In  section  6,  the  conclusions 
are  presented. 

2.  Dynamic  and  Thermodynamic  Integral  Equations 

Let  ( x ,  y )  and  z  be  the  horizontal  and  vertical  Cartesian  coordinates  with  the  constant  unit 
vectors  (ev,  ev)  and  ez.  The  circulation  model  is  established  on  the  base  of  the  two  approximations: 
hydrostatic  and  anelastic.  The  anelastic  approximation  is  to  assume  that  the  local  time  rate  of 
change  in  density  ( p  )  is  small;  the  continuity  equation  may  be  approximated  by  (Ogura  and 
Phillips,  1962) 

V  •  (pV)  =  o .  (1) 

Here,  V  =  ( u ,  v,  w)  is  the  velocity  vector  and  V  is  the  three-dimensional  gradient  operator.  The 
momentum  equation  is  given  by 

d(Py)  +  v  .  (pVV)  =  -Vp  +  V  •  (jiW)  +  F ,  (2a) 

dt 
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where  p  is  the  pressure,  //  is  the  eddy  viscosity.  F  represents  the  Coriolis  force  and  gravity  (body 
forces).  Let  f  be  a  scalar  representing  temperature,  salinity,  satisfying  the  advection-diffusion 


equation 


+  V  •  ( V  =  V  •  (zr^V  y/)  +  Fv, 


(2b) 


where  k  and  Fi/;  are  the  mixing  coefficient  and  the  forcing  term  for  /// .  The  eddy  viscosity  and 


vertical  mixing  coefficients  //  and  ^  using  the  level-2  turbulence  closure  (Mellor  and  Yamada, 
1974). 

Integration  of  (1)  on  a  finite  volume  £2  (Fig.  1)  leads  to 

I_qV  '(P'V)dQ  =  <j)  p\  •  nr/r , 
and  integration  of  (2a)  and  (2b)  leads  to 

•  ndT  =  -^  pndT  +  ij)  //VV  •  nr/r  +  j^FJQ , 


(3) 


+  =  ^ K  V y/ •  n<r/r  +  , 


(4a) 

(4b) 


where  T  is  the  boundary  of  Q  and  n  is  the  unit  vector  normal  to  T  (outward  positive).  The  two 
equations  (4a)  and  (4b)  are  quite  similar  in  tenns  of  using  the  finite  volume  method.  The  only 

difference  is  the  pressure  term,  -^pndT .  Equations  (4a)  and  (4b)  can  be  combined  into  one 


equation 


L  37^  +  * nclr  =  §r  ^  *  n< iT  +  JQ  F/ Q  ^ 


dt 


(4c) 


where  the  scalar  0is  one  of  ( pu,pv,if /)  and 


<p  =  pu,  K=ju,  AR  =  APx 
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(5) 


<P  =  pv,  *>=//,  AP0  =  APv, 

0  =  ^.  K+  =  K¥,  AP,=0, 

where  ( APx,AP  )  are  the  (x,  y)  components  of  the  pressure  gradient  force  on  the  volume  Q  . 

For  a  temporally  varying  finite  volume,  the  time  rate  of  change  of  the  volume  integrated  <p- 
value  is  computed  by 

|-(j)  <pdQ  =  ( £>  ^td&  +  </>^.  (6) 

dt  dt  dt 

Substitution  of  (6)  into  (4c)  leads  to 

^ \j>d&  -  </>^~  +  (f^V  •  ndT  =  j)rKf70»  ndT  +  +  A ,  (7a) 

which  is  the  basic  equation  for  (j)  (integral  0 -equation).  Time  integration  of  (7a)  from  t\  to  h  gives 

J  (f>{t2)dQ.  -^^(p{tx)dQ.  =  -Afij)  0(t*)X  •  ndT 

+Atd)  +  F,(f*)rfQ  +  AP,(f*)l  +  AfL— 1  ,  (7b) 

Jr  LJn  -I  y  dt  ),* 

where  A t  =  t2-tl,  and  t]  <t*<t2.  If  t*  =  t\,  the  scheme  is  explicit;  if  t*  =  h  ,  the  scheme  is  implicit. 

Adjustment  of  t*  may  lead  to  a  high-order  temporal  discrete  scheme.  The  last  term  in  the  right-hand 

side  of  (7b)  represents  the  temporal  volume  change  due  to  surface  elevation  fluctuation. 

3.  Discretization 

3.1.  Terrain-Following  Crystal  Grids 

Discretization  of  VOM  in  the  horizontal  directions 

x  =  x(i),  y  =  y(j)  ,  (8) 

is  similar  to  the  z-  and  <7  -coordinate  systems.  Discretization  of  VOM  in  the  vertical  direction. 


6 


z  =  z(i,j,  k,  t ), 


(9) 


varies  with  the  location  ( i,j ,  k )  and  time  (t).  Let  the  bottom  topography  is  represented  by 

z  =  -H(x,y), 


where  fl(x,  y )  is  assumed  to  be  a  singled  value  function  of  (x,  y).  Let  (ev,  ev)  be  the  constant  unit 
vectors  in  the  (x,  y)  directions. 

Four  types  of  finite  volumes  are  constructed.  Type-A  (Fig.  la)  has  four  vertical  lateral 
surfaces  perpendicular  to  er  or  ev, 


n  =  e 


nv  =  ev 


and  the  lower  surface  either  away  or  at  the  ocean  bottom.  The  four  vertices  (SW,  NW,  SE,  NE)  of 
the  upper  surface  are  away  from  the  bottom  topography.  The  upper  and  lower  slanted  surfaces  are 
trapezoids.  Type-B  has  one  vertex  of  the  upper  surface  at  the  bottom  topography  (say  vertex  SW). 
Since  H(x,  y)  is  a  single  valued  function,  it  is  all  land  below  the  vertex  SW.  The  lower  surface  of 
the  finite  volume  is  a  triangle  (no  longer  a  trapezoid).  Such  a  volume  (Fig.  lb)  still  has  four  vertical 
surfaces  and  two  slanted  surfaces  (both  triangles).  Type-C  has  two  vertices  of  the  upper  surface  at 
the  bottom  topography  (say  vertices  SW  and  NW).  It  is  all  land  below  the  vertices  SW  and  NW. 
The  lower  surface  shrinks  into  a  line.  The  finite  volume  has  two  slanted  surfaces  and  three  vertical 


surfaces  with  two  surfaces  perpendicular  to  nv  and  one  surface  perpendicular  to  n*  (Fig.  lc).  Type- 
D  has  three  vertices  of  the  upper  surface  at  the  bottom  topography  with  only  one  vertex  (say  NE)  in 
the  water.  The  vertices  NW  and  SE  are  at  the  bottom  boundary  and  the  vertex  SW  is  inside  the  land 
(not  shown  in  Fig.  Id).  The  lower  surface  shrinks  into  a  point.  The  finite  volume  has  two  slanted 
surfaces  and  two  vertical  surfaces  perpendicular  to  nx  and  nv  (Fig.  Id).  This  grid  system  is  called 
the  crystal  grid  due  to  the  crystal  shape  of  the  finite  volumes.  The  nonnal  unit  vector  ( n  s)  on  the 
slanted  surfaces  of  all  the  four  types  of  finite  volumes  is  calculated  by 
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n  = 


dz 


dz 


'  dx  n-r  dyny+*z 


\ 

L 

(  dz^ 

2 

(  3z^| 

— 

+ 

— 

+  1 . 


(ID 


The  mesh  characteristics  make  the  VOM  with  crystal  grid  (Fig.  1)  superior  to  the  z-  and  a  - 
coordinate  systems.  Fig.  2  shows  the  discretization  of  VOM  in  the  vertical  direction  differs  from 
both  z-  and  a  -  coordinate  models.  The  finite  volumes  with  the  crystal  grids  have  six  surfaces  for 
non-boundary  volumes  and  five  surfaces  for  boundary  volumes.  To  detennine  ns  in  a  finite  volume 
enclosed  by  (x„  x  ,+  /)  in  the  x-direction  and  by  (yj,  yj+i)  in  the  y-direction,  a  bi-linear  interpolation 
is  used  to  obtain  the  s-boundary  (i.e.,  slanted  boundary)  for  that  volume, 


z(x,y)  =  zn 


(X  ~  xM  )(y  -  yj+l )  (x- XM  )(y  -  y} ) 

- : - h  Z,  - — 


( V  - xM )(Vj  - yj+i )  12  ( xt  ~ xM )(Vj  -  Vj ) 

(x-x^iy-yj) 


+z  (x  ~  xi)(y  -  y  j+l) 

21  (xM  -  xtXyj  ~  yj+ 1)  22  (xi+l  -  v)(t/+i  -  yj) 


(12) 


where 

Zn  =  z(xi,yj),  z12  =  z(x„yJ+1),  z21  =  z(x,+1,yy),  z22  =  z(xM,yj+1)  . 

Substitution  of  (12)  into  (9)  determines  the  unit  vector  n  v  at  the  center  of  the  horizontal  cell,  (xj+1/2, 
Vj+1/2 )  =  [(*,-  +  xi+i  )/2,  (yj  +  yj+J  )/2], 

ns=(cos<9v,  cos(9v,  cos(9z),  (13) 

where 


cos  6 .  =  - 


zAxMn,yj+m) 


f 


,  2  .  2 
+  Zv+Zv 


cosO.=-ZyiXi+m,yj+m) 


1  +  K+Zy 


COSiT 


f 


I  2  .  2 

+  Z; v  +  Zy 


and  zv  =  dz  l  dx,  z  =  dz/dy  .  The  normal  gradient  of  any  variable  <j)  is  calculated  by 


!^  =  V**n,= 

dn. 


d(p  d(j)  d(j) 


dz  dx 


dy 


if 


,  2  .  2 
+  Z*+Zy 


(14) 
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3.2.  Flux  Computation 

The  surfaces  of  a  finite  volume  perpendicular  to  (nT,  nv,  nv)  are  called  the  x-,  y-  and  s- 
boundaries.  In  basic  equations  (4a)  and  (6),  the  fluxes  across  the  x-boundary  of  a  finite  volume  are 
computed  by 

/m=  f  .  .  py*nxdT  =  puAyAz,  fx  =  f  <pV  •  nxdY  =  tftuAyAz , 

J  x-boundary  T  J  x-boundary 

/,'=J\  ,  p;=pAyAz,  p;=0,  (is) 

J  x-boundary  T  T  rjy 

where  the  overbar  indicates  the  spatially  averaged  value  at  the  boundary  surface,  ( P*,P '*)  are  the 

two  components  of  the  horizontal  pressure  force  on  the  x-boundary.  The  fluxes  across  the  y- 
boundary  of  a  finite  volume  are  computed  by 


fm=  f  ,  ,  pV  •  ndT  =  pvAxAz,  //  =  [  </>Y  •  nvdT  =  ^vAxAz 

J  y-boundary  y  J  y-boundary  J 


y-boundary 


ft  =  f  ,  ,  •  n  vdT  =  ^-AxAz,  Py  =  pAxAz,  Pxy  =  0 . 

J  y-boundary  ■  T  Hi;  y 


J  y-boundary  v  '  y  r  y 


(16) 


where  ( P-  ,  P;'  )  are  the  two  components  of  the  horizontal  pressure  force  on  the  y-boundary.  The 


fluxes  across  the  s-boundary  of  a  finite  volume  are  computed  by 


/:=! 


pV  •  nsc/r  =  p(yv  —  u  z~~  —  v  ^-)AxAy  , 

s-boundary  Hy  n\? 


//={  ,  ,  0V»nsdr  =  0(w-u^-v^-)AxAy, 

T  J  s-boundary  dv  Hy 


d(j)  d(p  d(p  dz  d(pdz 

K,  — —dT  =  Ki—  — — - — — )AxAv  , 

s— boundary  0  0^  '  02  dx  0X  0  V  0V 


^  =  -pz^bxAy,  Py  =  ~P  %AxAy  ’ 


(17) 
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where  ( P/ ,  P/ )  are  the  two  components  of  the  horizontal  pressure  force  on  the  s-boundarv.  It  is 
noticed  that  the  derivatives  should  be  computed  by 


d(j)  _  8(j)  dz  8<p  d(j)  _  8<p  dz  8(j)  d(j)  _  8(j) 

dx  8x  dx  8z  dy  8y  dy  8z  dz  8z 


(18) 


where  ( 8(f)  /  8x,  8<j)  /  8 y,  80/  Sz  )  arc  computed  directly  from  the  grid  point  data. 

3.3.  Explicit  Scheme 

The  finite  volume  scheme  for  solving  continuity  equation  (3)  at  (i,j,  k )  is 


(19) 


The  finite  volume  scheme  for  solving  the  basic  equation  (7b)  is  explicit  if  t*  =  ti. 


vG  -  = m  -  /;  o' mu  fjaj +^,k) 


+ fl  (iJ,k-^)~  /;  (i,j,k  +  ^)-fTx(i-^,j,k)  +  fTx(i  +  ^,j,k)~  //  (i,j~,k) 

//  (i,j  +  i  ,k)  j,k  -  \)  +  f;  (i,  j,k  +  i)  +  F(i,  j,k)Q,.k  +  AP,  (i,  j,k)  +  At(<P 

222  ^  a? 


(20) 


where  n  is  the  time  step.  The  upper-most  finite  volumes  (k  =  1)  Q.  .j  changes  with  time  when  the 


surface  elevation  r/  varies, 


V  dt  Aj,i 


V  dt  Jtj  ’ 


The  finite  volumes  below  the  surface  Q.  . ,k  (k  >  1)  do  not  change  with  time. 


V  dt  )ULk 


0,  for  k  *  0 . 
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The  horizontal  pressure  gradient  force  ( AP  ,  APv )  on  the  finite  volume  is  computed  by 


a px = p;  a  ■ -  \ ,  j,  k )  -  p;  a + 1 ,  j,  k ) + p;  a,  j,  k-l-)-p:  a,  j,  k + , 


(21) 


(22) 


4.  Preliminary  Test 

4.1.  Test  Strategy 

Usually,  verification  of  a  new  numerical  model  should  be  divided  into  stages:  (1)  evaluating 
its  own  performance,  and  (2)  identifying  its  difference  from  the  existing  models.  Theoretically,  the 
performance  of  any  numerical  ocean  model  should  only  be  tested  against  analytical  or  known 
solutions.  For  coastal  ocean,  it  hardly  finds  any  analytical  solutions.  Known  solution  is  hard  to  find. 
Without  atmospheric  and  lateral  forcing,  the  ocean  that  is  initially  at  rest  should  be  at  rest  forever. 
Thus,  we  have  the  known  solutions  for  this  case, 

V  =  0,  |U0,  fU0.  (23) 

ox  oy 

The  seamount  test  case  (Beckmann  and  Haidvogel,  1993)  is  to  use  this  known  solution  for  model 
evaluation.  Any  nonzero  horizontal  pressure  gradients  (or  velocities)  obtained  from  the  numerical 
model  are  considered  errors. 

Several  advanced  test  cases  have  been  proposed  to  identify  the  model-model  difference  and 
the  sensitivity  to  the  choice  of  (say)  advection  algorithm,  such  as  gravitational  adjustment  of 
density  front,  residual  circulation  over  a  coastal  canyon,  combined  effects  of  topography  and 
stratification  (Haidvogel  and  Beckmann,  1998).  The  second  stage  test  should  be  conducted  after  the 
first  stage  test.  In  this  paper,  we  only  present  the  first  stage  evaluation  using  the  seamount  test  case. 


11 


In  regional  oceanic  (or  atmospheric)  prediction  models,  the  effects  of  bottom  topography 
must  be  taken  into  account  and  usually  the  terrain-following  sigma  coordinates  should  be  used  to 
imply  the  continuous  topography.  In  sigma  coordinates  the  water  column  is  divided  into  the  same 
number  of  grid  cells  regardless  of  the  depth.  Consider  2D  problems  for  mathematical  simplification. 
Let  ( x ,  z )  be  the  Cartesian  coordinates  and  ( x ,  a)  be  the  sigma  coordinates.  The  conventional 
relationships  between  z-  and  a  -coordinates  are  given  by 

x  =  x,  a=  Z  71  ,  (24) 

H  +  rj 

where  ij  is  the  surface  elevation.  Both  z  and  a  increase  vertically  upward  such  that  z  =  t) ,  <j  =  0  at 
the  surface  and  <7  =  -1 ,  z  =  -H  at  the  bottom.  The  horizontal  pressure  gradient  becomes  difference 
between  two  large  terms 

dJLj± - !_(CT^3)3p,  (25) 

dx  dx  H  +  a  dx  dx  d<7 

that  may  cause  large  truncation  error  at  steep  topography  (e.g.,  Gary  1973;  Haney  1991;  Mellor  et 
al.  1994;  McCalpin  1994;  Chu  and  Fan,  1997,  1998,  1999,  2000,  2001,  2003;  Song  1998).  Since 
the  horizontal  pressure  gradient  error  is  a  key  problem  in  the  terrain-following  ocean  models,  the 
first  step  of  the  VOM  test  should  be  the  evaluation  of  its  capability  on  reducing  the  horizontal 
pressure  gradient  error. 

4.2,  Seamount  Test  Case 

Suppose  a  seamount  located  inside  a  periodic  /-plane  (/  =  10"4  s'1)  channel  with  two  solid, 
free-slip  boundaries  along  constant  v.  Unforced  flow  over  seamount  in  the  presence  of  resting,  level 
isopycnals  is  an  idea  test  case  for  the  assessment  of  pressure  gradient  errors  in  simulating  stratified 
flow  over  topography.  The  flow  is  assumed  to  be  reentrant  (periodic)  in  the  along  channel 
coordinate  (i.e.,  x-axis).  The  seamount  test  case  is  chosen  to  test  the  performance  of  VOM. 
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The  domain  is  a  periodic  channel,  300  km  long  and  300  km  wide.  The  channel  has  a  far- 
field  depth  hmax  and  in  the  center  includes  an  isolated  Gaussian-shape  seamount  with  a  width  L 
and  an  amplitude  (Fig.  3), 


h(x,y)  =  h  li¬ 


ce  exp 


(x  —  Xq)2  +  ( v  —  v0)2 

1} 


(25) 


where  (x0,  v0)  are  the  longitude  and  latitude  of  the  seamount  center.  In  this  study,  we  use 

himx  =  4500  km,  L-  25  km,  a  -  0.9  .  (26) 

Suppose  that  the  background  fluid  is  at  rest  and  with  a  constant  salinity  (35  ppt)  and  an 
exponentially  stratified  initial  temperature, 

T(z)  =  5  +  1 5  exp(^— )  (unit:  °C )  ,  (27) 

where  Ht=  1000  m.  Since  the  fluid  is  initially  at  rest  and  the  density  field  is  independent  on  x  and  y, 
without  forcing  the  velocity  and  horizontal  pressure  gradient  should  be  zero.  Any  nonzero 
velocities  are  computational  errors. 


4.3.  Experiment  Setting 

A  a  -coordinate  finite  difference  model,  the  Princeton  Ocean  Model  (POM,  Blumberg  and 
Mellor,  1987),  is  implemented  for  the  seamount  test  case  using  horizontally  varying  grids  with  high 
resolution  over  the  seamount, 


(Ax),.  =  8  km 


1  -  0.5  sin 

'  in  \ 

UJJ 

(Ay)y  =  8  km 


(  •  ^ 

1  -  0.5  sin 

in 

M 

\  y ) 

j  =  1,2  ,...,Mv. 


(28) 
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where  Mx  =  My=  64.  The  VOM  with  the  same  physics,  parameterization,  and  horizontal  grids  as 
the  POM  is  also  implemented  for  seamount  test  case.  The  vertical  cross-sections  of  the  VOM  and 
POM  are  illustrated  in  Figs.  2b  and  2c.  The  time  steps  for  barotropic  and  baroclinic  modes  are  6  s 
and  180  s,  respectively. 

5.  Comparison  of  Finite  Difference  and  Finite  Volume  Schemes 

5.1.  Temporal  Variations  of  Error  Volume  Transport 

Both  cases  are  integrated  for  20  days  for  the  standard  test.  Fig.  4  displays  errors  in  the 
volume  transport  streamfunction  after  performing  20  days  of  integration  using  the  finite  volume 
and  finite  difference  schemes.  The  volume  transport  streamfunction  has  a  large-scale  eight-lobe 
pattern  centered  on  the  seamount.  The  errors  in  the  volume  transport  reduce  more  than  50%  from 
finite  difference  to  finite  volume  schemes.  On  the  20th  day,  the  errors  in  the  volume  transport 
varies  from  -56  to  84 x  10’3  Sv  using  the  POM  and  from  -28  to  45.5  x  10’3  Sv  using  the  VOM. 

5.2.  Temporal  Variations  of  Pressure  Gradient  Error 

Owing  to  a  very  large  number  of  calculations  performed,  we  discuss  the  results  exclusively 
in  terms  of  the  maximum  and  spatially  averaged  absolute  values  of  the  horizontal  pressure  gradient 
errors,  called  the  maximum  pressure  gradient  error  (PGmax)  and  the  mean  pressure  gradient  error 
(PGm).  Figure  5a  and  5b  show  the  time  evolution  of  PGmax  and  PGm  for  the  first  20  days  of 
integration  using  the  finite  difference  and  finite  volume  schemes.  Both  errors  increase  with  time, 
however,  they  are  10-15  times  smaller  using  the  finite  volume  scheme  than  using  the  finite 
difference  scheme.  For  example,  at  Day-10,  PGmax  =  33.42  x  10‘  N/m  using  the  finite  difference 
scheme  and  PGmax  =2.16x  10’9  N/m3  using  the  finite  volume  scheme;  PGm=  0.449  x  10'9  N/m3 
using  the  finite  difference  scheme  and  PGm  =  0.04  x  10’  N/m  using  the  finite  volume  scheme;  at 
Day-20,  PGmax  =58.41  x  10’9  N/m3  using  the  finite  difference  scheme  and  PGmax  =  4.18x  10’9  N/m3 
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using  the  finite  volume  scheme;  PGm  =1.596x  10'  N/m  using  the  finite  difference  scheme  and 
PGm  =  0.150x  10'  N/m  using  the  finite  volume  scheme. 

5.3.  Temporal  Variations  of  Error  Velocity 

Owing  to  a  very  large  number  of  calculations  perfonned,  we  discuss  the  results  exclusively 
in  tenns  of  the  maximum  and  spatially  averaged  absolute  values  of  the  spurious  velocity  generated 
by  the  pressure  gradient  errors,  called  the  peak  error  velocity  (Vp)  and  the  mean  error  velocity  (Vm). 
Figures  5c  and  5d  show  the  time  evolution  of  the  mean  and  peak  error  velocity  for  the  first  20  days 
of  integration  using  the  finite  difference  and  finite  volume  schemes.  Both  peak  and  mean  error 
velocities  increase  with  time,  however,  they  are  4  times  smaller  using  the  finite  volume  scheme 
than  using  the  finite  difference  scheme.  For  example,  at  Day- 10,  Vp  =  2.61  cm/s  using  the  finite 
difference  scheme  and  Vp  =0.57  cm/s  using  the  finite  volume  scheme;  Vm  =  0.054  cm/s  using  the 
finite  difference  scheme  and  Vm  =0.015  cm/s  using  the  finite  volume  scheme;  at  Day-20,  Vp=  4.25 
cm/s  using  the  finite  difference  scheme  and  Vp  =  0.98  cm/s  using  the  finite  volume  scheme;  Vm  = 
1.033  cm/s  using  the  finite  difference  scheme  and  Vm  =  0.028  cm/s  using  the  finite  volume  scheme. 

6.  Conclusions 

(1)  A  three-dimensional,  finite  volume  ocean  circulation  model  has  been  developed.  The 
basic  equations  are  transformed  from  differential  into  integral  forms  using  the  hydrostatic  and 
anelastic  approximations  and  solved  for  finite  volumes  (rather  than  grid  points)  with  the  flux 
conservation  enforced  on  arbitrarily  meshes.  This  model  has  great  flexibility  in  establishing  model 
grids. 

(2)  A  terrain-following  grid  discretization  is  proposed  such  that  the  vertical  surfaces  are 
perpendicular  to  horizontal  (x,  y)  axes;  and  the  slanted  surfaces  follow  the  ocean  surface  and 
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bottom.  Such  a  grid  system  reveals  a  superior  feature  than  both  z-  and  a  -coordinate  systems  over 
abrupt  bottom  topography. 

(3)  Seamount  test  case  is  the  first  step  to  show  the  value-added  of  using  finite  volume 
scheme.  The  second-order  finite  volume  scheme  leads  to  a  drastic  error  reduction  comparing  to  the 
second-order  finite  difference  scheme  using  POM. 

(4)  It  is  noticed  that  the  seamount  test  case  presented  here  is  preliminary.  More  cases  should 
be  conducted  in  the  future  for  testing  the  difference  between  VOM  and  the  existing  ocean  models 
using  gravitational  adjustment  of  a  density  front,  residual  circulation  over  a  coastal  canyon,  and 
combined  effects  of  topography  and  stratification. 
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Figure  Legends 


Fig.  1.  Four  types  of  finite  volumes:  (a)  Type-A  with  four  vertices  (SW,  NW,  SE,  NE)  of  the  upper 
surface  away  from  the  bottom  topography,  (b)  Type-B  with  one  vertex  (SW)  of  the  upper  surface 
at  the  bottom  topography,  (c)  Type-C  with  two  vertices  (SW  and  NW)  of  the  upper  surface  at  the 
bottom  topography,  and  (d)  Type-D  with  three  vertices  of  the  upper  surface  at  the  bottom 
topography  with  only  one  vertex  (NE)  in  the  water. 

Fig.  2.  Comparison  among  (a)  z-coordinate,  (b)  finite  volume,  and  (c)  a  -coordinate  systems. 

Fig.  3.  Seamount  geometry. 

Fig.  4.  Volume  transport  streamfunction  (Sv)  at  day-20  using  the  finite  difference  and  finite  volume 
schemes. 

Fig.  5.  Comparison  between  the  finite  difference  and  finite  volume  schemes  on  temporal  variations 
of  (a)  maximum  pressure  gradient  error  (N/m  ),  (b)  mean  pressure  gradient  error  (N/m  ),  (c)  peak 
error  velocity  (m/s),  and  (d)  mean  error  velocity  (m/s). 
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Fig.  1.  Four  types  of  finite  volumes:  (a)  Type-A  with  four  vertices  (SW,  NW,  SE,  NE)  of  the  upper 
surface  away  from  the  bottom  topography,  (b)  Type-B  with  one  vertex  (SW)  of  the  upper  surface 
at  the  bottom  topography,  (c)  Type-C  with  two  vertices  (SW  and  NW)  of  the  upper  surface  at  the 
bottom  topography,  and  (d)  Type-D  with  three  vertices  of  the  upper  surface  at  the  bottom 
topography  with  only  one  vertex  (NE)  in  the  water. 
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Fig.  2.  Comparison  among  (a)  z-coordinate,  (b)  finite  volume,  and  (c)  o  -coordinate  systems. 
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Fig.  3.  Seamount  geometry. 
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Fig.  4.  Volume  transport  streamfunction  (Sv)  at  day-20  using  the  finite  difference  and  finite  volume 
schemes. 
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Fig.  5.  Comparison  between  the  finite  difference  and  finite  volume  schemes  on  temporal  variations 
of  (a)  maximum  pressure  gradient  error  (N/m3),  (b)  mean  pressure  gradient  error  (N/m3),  (c)  peak 
error  velocity  (m/s),  and  (d)  mean  error  velocity  (m/s). 


24 


25 


